## Base RF model
set.seed(42)
library(ggplot2) # for autoplot() generic
library(dplyr)
library(sf)
library(stars)
library(tmap)
dat_sf <- st_read("./data/glacier_clim.shp")
## Reading layer `glacier_clim' from data source `/Users/u0784726/Dropbox/Data/devtools/glacier_mb/data/glacier_clim.shp' using driver `ESRI Shapefile'
## Simple feature collection with 95086 features and 31 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 67.47845 ymin: 27.4918 xmax: 103.8915 ymax: 45.35089
## Geodetic CRS: GCS_unknown
Load and crop basemap
basemap <- read_stars("~/Dropbox/Data/summer/natural_earth/HYP_50M_SR_W/HYP_50M_SR_W.tif")
basemap2 <- st_crop(basemap, st_bbox(dat_sf))
basemap2 <- st_rgb(basemap2)
t2mt2m <- read_stars("~/Dropbox/Data/summer/era5/t2m_2001_2020_amjjas.nc")
st_crs(t2m) <- 4326
t2m <- st_crop(t2m, st_bbox(dat_sf))
t2m <- st_as_stars(t2m)
Make average 2006-2010
t2m_2008 <- t2m %>%
slice("time", 6:10) %>%
st_apply(MARGIN = 1:2, FUN = mean)
Make average 2016-2020
t2m_2018 <- t2m %>%
slice("time", 16:20) %>%
st_apply(MARGIN = 1:2, FUN = mean)
Make delta
t2m_d <- t2m_2018 - t2m_2008
tptp <- read_stars("~/Dropbox/Data/summer/era5/tp_2001_2020_ann.nc")
st_crs(tp) <- 4326
tp <- st_crop(tp, st_bbox(dat_sf))
tp <- st_as_stars(tp)
Make average 2006-2010
tp_2008 <- tp %>%
slice("time", 6:10) %>%
st_apply(MARGIN = 1:2, FUN = mean)
Make average 2016-2020
tp_2018 <- tp %>%
slice("time", 16:20) %>%
st_apply(MARGIN = 1:2, FUN = mean)
Make delta
tp_d <- tp_2018 - tp_2008
t2mx <- c(t2m_2008, t2m_2018)
names(x) <- c("t2m_2008", "t2m_2018")
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(x) +
tm_raster(palette = "-inferno",
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\)
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(t2m_d) +
tm_raster(alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\) with mass balance
m1 <- tm_shape(t2m_d) +
tm_raster(alpha = 0.75,
style = "quantile") +
tm_shape(dat_sf) +
tm_symbols(col = "mb_mwea",
size = 0.25,
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
tpx <- c(tp_2008, tp_2018)
names(x) <- c("tp_2008", "tp_2018")
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(x) +
tm_raster(palette = "YlGnBu",
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\)
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(tp_d) +
tm_raster(alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\) with mass balance
m1 <- tm_shape(tp_d) +
tm_raster(alpha = 0.75,
style = "quantile") +
tm_shape(dat_sf) +
tm_symbols(col = "mb_mwea",
size = 0.25,
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)